! -----------------------------------------------------------------------
! File:             Basic_matrix.f90
! Purpose:          Definition of basic vectorized operators
! Contains:
! Revision History: Ver. 1.0 May. 2007 Alberto Scotti
! -----------------------------------------------------------------------

MODULE Basic_matrix
USE ElliDef, ONLY: Precision
IMPLICIT NONE
SAVE 
PRIVATE

PUBLIC gc2e,ge2c,e2c,c2e,gc2eb,c2eb,strip,zeroes,zeroi,dress


REAL (KIND=Precision), PARAMETER :: zero = 0.0D0,half=0.5D0,one=1.0d0
CONTAINS

SUBROUTINE gc2e(a,ja,ia,nNonZeros,delta,dir,nx,ny,nz)
! Defines and fills in the sparse matrix corresponding to a gradient 
! in the direction specified by dir, applied to a cell-centered field 
! and returned on cell edges (c2e) along which the grid spacing is
! delta. The sparse matrix is (nx+I)*(ny+J)*(nz+K) X nx*ny*nz where
! I=1 iff dir=1, J=1 iff dir=2 and K=1 iff dir=3 
! NOTE THAT DR BOUNDARY CONDITIONS ARE ASSUMED. IT IS UP TO THE USER
! IN THE CALLING ROUTINE TO MAKE SURE THAT BOUNDARY VALUES ARE MODIFIED
! ACCORDINGLY. NOTE THAT FOR Z DERIVATIVES, THE LOWER BOUNDARY CONTAINS
! PLACE HOLDERS TO ACCOMODATE COMPLEX BOUNDARIES.
! Declarations
INTEGER, INTENT (IN) :: nx,ny,nz,dir
REAL (KIND=Precision), INTENT (IN) :: delta
INTEGER, INTENT (OUT) :: nNonZeros
INTEGER, POINTER, DIMENSION (:) :: ia,ja
REAL (KIND=Precision), POINTER, DIMENSION (:) :: a

INTEGER :: i,j,k,l,m,d(3,3)


! release memory
IF(ASSOCIATED(a)) THEN
   DEALLOCATE(a); NULLIFY(a)
ENDIF
IF(ASSOCIATED(ia)) THEN
   DEALLOCATE(ia); NULLIFY(ia)
ENDIF
IF(ASSOCIATED(ja)) THEN
   DEALLOCATE(ja); NULLIFY(ja)
ENDIF


SELECT CASE (dir)

CASE (1) ! gradient is along x direction
! determine size of matrix and allocate space
nNonZeros=nz*ny*((nx-1)*2+2)
ALLOCATE(a(nNonZeros),ja(nNonZeros),ia((nx+1)*ny*nz+1))

nNonZeros=1
   DO k=1,nz
      DO j=1,ny
         DO i=1,nx+1
            l=(k-1)*ny*nx+(j-1)*nx+i ! vectorized column coordinate
            m=(k-1)*ny*(nx+1)+(j-1)*(nx+1)+i ! vector. row coordinate
            IF(i==1) THEN
               a(nNonZeros:nNonZeros)=1./delta
               ja(nNonzeros:nNonZeros)= l
               ia(m)=nNonZeros
               nNonZeros=nNonZeros+1 ! reposition pointer to next row
            ELSEIF(i==nx+1) THEN
               a(nNonZeros:nNonZeros)=-1./delta
               ja(nNonzeros:nNonZeros)= l-1
               ia(m)=nNonZeros
               nNonZeros=nNonZeros+1 ! reposition pointer to next row
            ELSE
               a(nNonZeros:nNonZeros+1)=(/-1./delta,1./delta/)
               ja(nNonZeros:nNonZeros+1)=(/l-1,l/)
               ia(m)=nNonzeros
               nNonZeros=nNonZeros+1+1
            ENDIF
         ENDDO
      ENDDO
   ENDDO
nNonZeros=nNonZeros-1;ia(m+1)=nNonZeros+1
CASE (2) ! gradient is along y direction
nNonZeros=nz*nx*((ny-1)*2+2)
ALLOCATE(a(nNonZeros),ja(nNonZeros),ia(nx*(ny+1)*nz+1))

nNonZeros=1
   DO k=1,nz
      DO j=1,ny+1
         DO i=1,nx
            l=(k-1)*ny*nx+(j-1)*nx+i ! vectorized coordinate 
            m=(k-1)*(ny+1)*nx+(j-1)*nx+i ! vector. row coordinate
            IF(j==1) THEN
               a(nNonZeros:nNonZeros)=1./delta
               ja(nNonzeros:nNonZeros)= l
               ia(m)=nNonZeros
               nNonZeros=nNonZeros+1 ! reposition pointer to next row
            ELSEIF(j==ny+1) THEN
               a(nNonZeros:nNonZeros)=-1./delta
               ja(nNonzeros:nNonZeros)= l-nx
               ia(m)=nNonZeros
               nNonZeros=nNonZeros+1
            ELSE
               a(nNonZeros:nNonZeros+1)=(/-1./delta,1./delta/)
               ja(nNonZeros:nNonZeros+1)=(/l-nx,l/)
               ia(m)=nNonzeros
               nNonZeros=nNonZeros+1+1
            ENDIF
         ENDDO
      ENDDO
   ENDDO
nNonZeros=nNonZeros-1;ia(m+1)=nNonZeros+1

CASE (3) ! gradient is along z direction
   nNonZeros=ny*nx*((nz-1)*2+1)+(nx-2)*(ny-2)*5+(nx-2)*8+(ny-2)*8+12
   ALLOCATE(a(nNonZeros),ja(nNonZeros),ia(nx*ny*(nz+1)+1))
   
   nNonZeros=1
   DO k=1,nz+1
      DO j=1,ny
         DO i=1,nx
            l=(k-1)*ny*nx+(j-1)*nx+i ! vectorized coordinate 
            m=(k-1)*ny*nx+(j-1)*nx+i ! vector. row coordinate
            IF(k==1) THEN ! at bottom, we add space holders for non-diagonal
                          ! metric induced boundary conditions
               IF(MOD(i-1,nx-1)+MOD(j-1,ny-1)==0) THEN ! vertices
                  IF(i==1.AND.j==1) THEN
                     a(nNonZeros:nNonZeros+2)=(/1./delta,zero,zero/)
                     ja(nNonZeros:nNonZeros+2)=(/l,l+1,l+nx/)
                     ia(m)=nNonzeros
                     nNonZeros=nNonZeros+2+1
                  ELSEIF(i==1.AND.j==ny) THEN
                     a(nNonZeros:nNonZeros+2)=(/zero,1./delta,zero/)
                     ja(nNonZeros:nNonZeros+2)=(/l-nx,l,l+1/)
                     ia(m)=nNonzeros
                     nNonZeros=nNonZeros+2+1
                  ELSEIF(i==nx.AND.j==1) THEN
                     a(nNonZeros:nNonZeros+2)=(/zero,1./delta,zero/)
                     ja(nNonZeros:nNonZeros+2)=(/l-1,l,l+nx/)
                     ia(m)=nNonzeros
                     nNonZeros=nNonZeros+2+1
                  ELSE
                      a(nNonZeros:nNonZeros+2)=(/zero,zero,1./delta/)
                     ja(nNonZeros:nNonZeros+2)=(/l-nx,l-1,l/)
                     ia(m)=nNonzeros
                     nNonZeros=nNonZeros+2+1
                  ENDIF
               ELSEIF(MOD(i-1,nx-1)==0) THEN
                  IF(i==1) THEN
                     a(nNonZeros:nNonZeros+3)=(/zero,1./delta,zero,zero/)
                     ja(nNonZeros:nNonZeros+3)=(/l-nx,l,l+1,l+nx/)
                     ia(m)=nNonzeros
                     nNonZeros=nNonZeros+3+1
                  ELSE
                     a(nNonZeros:nNonZeros+3)=(/zero,zero,1./delta,zero/)
                     ja(nNonZeros:nNonZeros+3)=(/l-nx,l-1,l,l+nx/)
                     ia(m)=nNonzeros
                     nNonZeros=nNonZeros+3+1
                  ENDIF
               ELSEIF(MOD(j-1,ny-1)==0) THEN
                  IF(j==1) THEN
                     a(nNonZeros:nNonZeros+3)=(/zero,1./delta,zero,zero/)
                     ja(nNonZeros:nNonZeros+3)=(/l-1,l,l+1,l+nx/)
                     ia(m)=nNonzeros
                     nNonZeros=nNonZeros+3+1
                  ELSE
                     a(nNonZeros:nNonZeros+3)=(/zero,zero,1./delta,zero/)
                     ja(nNonZeros:nNonZeros+3)=(/l-nx,l-1,l,l+1/)
                     ia(m)=nNonzeros
                     nNonZeros=nNonZeros+3+1
                  ENDIF
               ELSE
                  a(nNonZeros:nNonZeros+4)=(/zero,zero,1./delta,zero,zero/)
                  ja(nNonzeros:nNonZeros+4)= (/l-nx,l-1,l,l+1,l+nx/)
                  ia(m)=nNonZeros
                  nNonZeros=nNonZeros+4+1 ! reposition pointer to next row
               ENDIF
            ELSEIF(k==nz+1) THEN
               a(nNonZeros:nNonZeros)=-1./delta
               ja(nNonzeros:nNonZeros)= l-nx*ny
               ia(m)=nNonZeros
               nNonZeros=nNonZeros+1
            ELSE
               a(nNonZeros:nNonZeros+1)=(/-1./delta,1./delta/)
               ja(nNonZeros:nNonZeros+1)=(/l-nx*ny,l/)
               ia(m)=nNonzeros
               nNonZeros=nNonZeros+1+1
            ENDIF
         ENDDO
      ENDDO
   ENDDO
nNonZeros=nNonZeros-1;ia(m+1)=nNonZeros+1

END SELECT

IF(size(a)/=nNonzeros) THEN
   PRINT*,'nNonzeros is not correct in gc2e'
   STOP
ENDIF
IF(MAXVAL(ja)>nx*ny*nz) THEN
   PRINT*,'problems with definition of matrix in gc2e'
   PRINT*,'columns exceed size'
   STOP
END IF

END SUBROUTINE gc2e

SUBROUTINE ge2c(a,ja,ia,nNonZeros,delta,dir,nx,ny,nz)
! Defines and fills in the sparse matrix corresponding to a gradient 
! in the direction specified by dir, applied on a edge field and returned 
! on a cell centered field along which the grid spacing is
! delta. The sparse matrix is nx*ny*nz X (nx+I)*(ny+J)*(nz+K) where
! I=1 iff dir=1, J=1 iff dir=2 and K=1 iff dir=3 
! Declarations
INTEGER, INTENT (IN) :: nx,ny,nz,dir
REAL (KIND=Precision), INTENT (IN) :: delta
INTEGER, INTENT (OUT) :: nNonZeros
INTEGER, POINTER, DIMENSION (:) :: ia,ja
REAL (KIND=Precision), POINTER, DIMENSION (:) :: a

INTEGER :: i,j,k,l,m,ncol


! release memory
IF(ASSOCIATED(a)) THEN
   DEALLOCATE(a); NULLIFY(a)
ENDIF
IF(ASSOCIATED(ia)) THEN
   DEALLOCATE(ia); NULLIFY(ia)
ENDIF
IF(ASSOCIATED(ja)) THEN
   DEALLOCATE(ja); NULLIFY(ja)
ENDIF


SELECT CASE (dir)

CASE (1) ! gradient is along x direction
! determine size of matrix and allocate space
nNonZeros=nz*ny*nx*2;ncol=(nx+1)*ny*nz
ALLOCATE(a(nNonZeros),ja(nNonZeros),ia((nx)*ny*nz+1))

nNonZeros=1
   DO k=1,nz
      DO j=1,ny
         DO i=1,nx

            l=(k-1)*ny*(nx+1)+(j-1)*(nx+1)+i ! vectorized coordinate column 
            m=(k-1)*ny*nx+(j-1)*nx+i ! row
               a(nNonZeros:nNonZeros+1)=(/-1./delta,1./delta/)
               ja(nNonZeros:nNonZeros+1)=(/l,l+1/)
               ia(m)=nNonzeros
               nNonZeros=nNonZeros+1+1

            
         ENDDO
      ENDDO
   ENDDO
nNonZeros=nNonZeros-1;ia(m+1)=nNonZeros+1
CASE (2) ! gradient is along y direction
nNonZeros=nz*nx*ny*2;ncol=nx*(ny+1)*nz
ALLOCATE(a(nNonZeros),ja(nNonZeros),ia(nx*(ny)*nz+1))

nNonZeros=1
   DO k=1,nz
      DO j=1,ny
         DO i=1,nx

            l=(k-1)*(ny+1)*(nx)+(j-1)*(nx)+i ! vectorized coordinate column 
            m=(k-1)*ny*nx+(j-1)*nx+i ! row
               a(nNonZeros:nNonZeros+1)=(/-1./delta,1./delta/)
               ja(nNonZeros:nNonZeros+1)=(/l,l+nx/)
               ia(m)=nNonzeros
               nNonZeros=nNonZeros+1+1

         ENDDO
      ENDDO
   ENDDO
nNonZeros=nNonZeros-1;ia(m+1)=nNonZeros+1

CASE (3) ! gradient is along z direction
nNonZeros=nz*nx*ny*2;ncol=nx*ny*(nz+1)
ALLOCATE(a(nNonZeros),ja(nNonZeros),ia(nx*(ny)*nz+1))

nNonZeros=1
   DO k=1,nz
      DO j=1,ny
         DO i=1,nx

            l=(k-1)*ny*(nx)+(j-1)*(nx)+i ! vectorized coordinate column 
            m=(k-1)*ny*nx+(j-1)*nx+i ! row
               a(nNonZeros:nNonZeros+1)=(/-1./delta,1./delta/)
               ja(nNonZeros:nNonZeros+1)=(/l,l+ny*nx/)
               ia(m)=nNonzeros
               nNonZeros=nNonZeros+1+1

         ENDDO
      ENDDO
   ENDDO
nNonZeros=nNonZeros-1;ia(m+1)=nNonZeros+1

END SELECT

IF(size(a)/=nNonzeros) THEN
   PRINT*,'nNonzeros is not correct in ge2c'
   PRINT*,'size of a is',size(a)
   PRINT*,'nNonZeros =',nNonZeros
   STOP
ENDIF
IF(MAXVAL(ja)>ncol) THEN
   PRINT*,'problems with definition of matrix in ge2c'
   PRINT*,'columns exceed size'
   STOP
END IF


END SUBROUTINE ge2c

SUBROUTINE e2c(a,ja,ia,nNonZeros,dir,nx,ny,nz)
! Defines and fills in the sparse matrix corresponding to a 
! edge-to-cell operation 
! in the direction specified by dir, applied on a edge field and returned 
! on a cell centered field.
! The sparse matrix is nx*ny*nz X (nx+I)*(ny+J)*(nz+K) where
! I=1 iff dir=1, J=1 iff dir=2 and K=1 iff dir=3 
! Declarations
INTEGER, INTENT (IN) :: nx,ny,nz,dir
INTEGER, INTENT (OUT) :: nNonZeros
INTEGER, POINTER, DIMENSION (:) :: ia,ja
REAL (KIND=Precision), POINTER, DIMENSION (:) :: a

INTEGER :: i,j,k,l,m,ncol


! release memory
IF(ASSOCIATED(a)) THEN
   DEALLOCATE(a); NULLIFY(a)
ENDIF
IF(ASSOCIATED(ia)) THEN
   DEALLOCATE(ia); NULLIFY(ia)
ENDIF
IF(ASSOCIATED(ja)) THEN
   DEALLOCATE(ja); NULLIFY(ja)
ENDIF


SELECT CASE (dir)

CASE (1) ! average is along x direction
! determine size of matrix and allocate space
nNonZeros=nz*ny*nx*2;ncol=(nx+1)*ny*nz
ALLOCATE(a(nNonZeros),ja(nNonZeros),ia((nx)*ny*nz+1))

nNonZeros=1
   DO k=1,nz
      DO j=1,ny
         DO i=1,nx

            l=(k-1)*ny*(nx+1)+(j-1)*(nx+1)+i ! vectorized coordinate column 
            m=(k-1)*ny*nx+(j-1)*nx+i ! row
               a(nNonZeros:nNonZeros+1)=(/half,half/)
               ja(nNonZeros:nNonZeros+1)=(/l,l+1/)
               ia(m)=nNonzeros
               nNonZeros=nNonZeros+1+1

            
         ENDDO
      ENDDO
   ENDDO
nNonZeros=nNonZeros-1;ia(m+1)=nNonZeros+1
CASE (2) ! gradient is along y direction
nNonZeros=nz*nx*ny*2;ncol=nx*(ny+1)*nz
ALLOCATE(a(nNonZeros),ja(nNonZeros),ia(nx*(ny)*nz+1))

nNonZeros=1
   DO k=1,nz
      DO j=1,ny
         DO i=1,nx

            l=(k-1)*(ny+1)*(nx)+(j-1)*(nx)+i ! vectorized coordinate column 
            m=(k-1)*ny*nx+(j-1)*nx+i ! row
               a(nNonZeros:nNonZeros+1)=(/half,half/)
               ja(nNonZeros:nNonZeros+1)=(/l,l+nx/)
               ia(m)=nNonzeros
               nNonZeros=nNonZeros+1+1

         ENDDO
      ENDDO
   ENDDO
nNonZeros=nNonZeros-1;ia(m+1)=nNonZeros+1

CASE (3) ! gradient is along z direction
nNonZeros=nz*nx*ny*2;ncol=nx*ny*(nz+1)
ALLOCATE(a(nNonZeros),ja(nNonZeros),ia(nx*(ny)*nz+1))

nNonZeros=1
   DO k=1,nz
      DO j=1,ny
         DO i=1,nx

            l=(k-1)*ny*(nx)+(j-1)*(nx)+i ! vectorized coordinate column 
            m=(k-1)*ny*nx+(j-1)*nx+i ! row
               a(nNonZeros:nNonZeros+1)=(/half,half/)
               ja(nNonZeros:nNonZeros+1)=(/l,l+ny*nx/)
               ia(m)=nNonzeros
               nNonZeros=nNonZeros+1+1

         ENDDO
      ENDDO
   ENDDO
nNonZeros=nNonZeros-1;ia(m+1)=nNonZeros+1

END SELECT

IF(size(a)/=nNonzeros) THEN
   PRINT*,'nNonzeros is not correct in e2c'
   STOP
ENDIF
IF(MAXVAL(ja)>ncol) THEN
   PRINT*,'problems with definition of matrix in e2c'
   PRINT*,'columns exceed size'
   STOP
END IF


END SUBROUTINE e2c

SUBROUTINE c2e(a,ja,ia,nNonZeros,dir,nx,ny,nz)
! Defines and fills in the sparse matrix corresponding to an averaging 
! in the direction specified by dir, applied to a cell-centered field 
! and returned on cell edges (c2e) along which the grid spacing is
! delta. The sparse matrix is (nx+I)*(ny+J)*(nz+K) X nx*ny*nz where
! I=1 iff dir=1, J=1 iff dir=2 and K=1 iff dir=3 
! 
! Declarations
INTEGER, INTENT (IN) :: nx,ny,nz,dir
INTEGER, INTENT (OUT) :: nNonZeros
INTEGER, POINTER, DIMENSION (:) :: ia,ja
REAL (KIND=Precision), POINTER, DIMENSION (:) :: a

INTEGER :: i,j,k,l,m


! release memory
IF(ASSOCIATED(a)) THEN
   DEALLOCATE(a); NULLIFY(a)
ENDIF
IF(ASSOCIATED(ia)) THEN
   DEALLOCATE(ia); NULLIFY(ia)
ENDIF
IF(ASSOCIATED(ja)) THEN
   DEALLOCATE(ja); NULLIFY(ja)
ENDIF


SELECT CASE (dir)

CASE (1) ! average is along x direction
! determine size of matrix and allocate space
nNonZeros=nz*ny*((nx-1)*2+2)
ALLOCATE(a(nNonZeros),ja(nNonZeros),ia((nx+1)*ny*nz+1))

nNonZeros=1
   DO k=1,nz
      DO j=1,ny
         DO i=1,nx+1
            l=(k-1)*ny*nx+(j-1)*nx+i ! vectorized coordinate 
            m=(k-1)*ny*(nx+1)+(j-1)*(nx+1)+i ! vector. row index
            IF(i==1) THEN
               a(nNonZeros:nNonZeros)=half
               ja(nNonzeros:nNonZeros)= l
               ia(m)=nNonZeros
               nNonZeros=nNonZeros+1 ! reposition pointer to next row
            ELSEIF(i==nx+1) THEN
               a(nNonZeros:nNonZeros)=half
               ja(nNonzeros:nNonZeros)= l-1
               ia(m)=nNonZeros
               nNonZeros=nNonZeros+1
            ELSE
               a(nNonZeros:nNonZeros+1)=(/half,half/)
               ja(nNonZeros:nNonZeros+1)=(/l-1,l/)
               ia(m)=nNonzeros
               nNonZeros=nNonZeros+1+1
            ENDIF
         ENDDO
      ENDDO
   ENDDO
nNonZeros=nNonZeros-1;ia(m+1)=nNonZeros+1
CASE (2) ! gradient is along y direction
nNonZeros=nz*nx*((ny-1)*2+2)
ALLOCATE(a(nNonZeros),ja(nNonZeros),ia(nx*(ny+1)*nz+1))

nNonZeros=1
   DO k=1,nz
      DO j=1,ny+1
         DO i=1,nx
            l=(k-1)*ny*nx+(j-1)*nx+i ! vectorized coordinate 
            m=(k-1)*(ny+1)*(nx)+(j-1)*(nx)+i ! vector. row index
            IF(j==1) THEN
               a(nNonZeros:nNonZeros)=half
               ja(nNonzeros:nNonZeros)= l
               ia(m)=nNonZeros
               nNonZeros=nNonZeros+1 ! reposition pointer to next row
            ELSEIF(j==ny+1) THEN
               a(nNonZeros:nNonZeros)=half
               ja(nNonzeros:nNonZeros)= l-nx
               ia(m)=nNonZeros
               nNonZeros=nNonZeros+1
            ELSE
               a(nNonZeros:nNonZeros+1)=(/half,half/)
               ja(nNonZeros:nNonZeros+1)=(/l-nx,l/)
               ia(m)=nNonzeros
               nNonZeros=nNonZeros+1+1
            ENDIF
         ENDDO
      ENDDO
   ENDDO
nNonZeros=nNonZeros-1;ia(m+1)=nNonZeros+1

CASE (3) ! gradient is along z direction
   nNonZeros=ny*nx*((nz-1)*2+1)+(nx-2)*(ny-2)*5+(nx-2)*8+(ny-2)*8+12
ALLOCATE(a(nNonZeros),ja(nNonZeros),ia(nx*ny*(nz+1)+1))

nNonZeros=1
   DO k=1,nz+1
      DO j=1,ny
         DO i=1,nx
            l=(k-1)*ny*nx+(j-1)*nx+i ! vectorized coordinate 
            m=(k-1)*ny*nx+(j-1)*nx+i ! vector. row index
            IF(k==1) THEN ! at bottom, we add space holders for non-diagonal
                          ! metric induced boundary conditions
               IF(MOD(i-1,nx-1)+MOD(j-1,ny-1)==0) THEN ! vertices
                  IF(i==1.AND.j==1) THEN
                     a(nNonZeros:nNonZeros+2)=(/half,zero,zero/)
                     ja(nNonZeros:nNonZeros+2)=(/l,l+1,l+nx/)
                     ia(m)=nNonzeros
                     nNonZeros=nNonZeros+2+1
                  ELSEIF(i==1.AND.j==ny) THEN
                     a(nNonZeros:nNonZeros+2)=(/zero,half,zero/)
                     ja(nNonZeros:nNonZeros+2)=(/l-nx,l,l+1/)
                     ia(m)=nNonzeros
                     nNonZeros=nNonZeros+2+1
                  ELSEIF(i==nx.AND.j==1) THEN
                     a(nNonZeros:nNonZeros+2)=(/zero,half,zero/)
                     ja(nNonZeros:nNonZeros+2)=(/l-1,l,l+nx/)
                     ia(m)=nNonzeros
                     nNonZeros=nNonZeros+2+1
                  ELSE
                      a(nNonZeros:nNonZeros+2)=(/zero,zero,half/)
                     ja(nNonZeros:nNonZeros+2)=(/l-nx,l-1,l/)
                     ia(m)=nNonzeros
                     nNonZeros=nNonZeros+2+1
                  ENDIF
               ELSEIF(MOD(i-1,nx-1)==0) THEN
                  IF(i==1) THEN
                     a(nNonZeros:nNonZeros+3)=(/zero,half,zero,zero/)
                     ja(nNonZeros:nNonZeros+3)=(/l-nx,l,l+1,l+nx/)
                     ia(m)=nNonzeros
                     nNonZeros=nNonZeros+3+1
                  ELSE
                     a(nNonZeros:nNonZeros+3)=(/zero,zero,half,zero/)
                     ja(nNonZeros:nNonZeros+3)=(/l-nx,l-1,l,l+nx/)
                     ia(m)=nNonzeros
                     nNonZeros=nNonZeros+3+1
                  ENDIF
               ELSEIF(MOD(j-1,ny-1)==0) THEN
                  IF(j==1) THEN
                     a(nNonZeros:nNonZeros+3)=(/zero,half,zero,zero/)
                     ja(nNonZeros:nNonZeros+3)=(/l-1,l,l+1,l+nx/)
                     ia(m)=nNonzeros
                     nNonZeros=nNonZeros+3+1
                  ELSE
                     a(nNonZeros:nNonZeros+3)=(/zero,zero,half,zero/)
                     ja(nNonZeros:nNonZeros+3)=(/l-nx,l-1,l,l+1/)
                     ia(m)=nNonzeros
                     nNonZeros=nNonZeros+3+1
                  ENDIF
               ELSE
                  a(nNonZeros:nNonZeros+4)=(/zero,zero,half,zero,zero/)
                  ja(nNonzeros:nNonZeros+4)= (/l-nx,l-1,l,l+1,l+nx/)
                  ia(m)=nNonZeros
                  nNonZeros=nNonZeros+4+1 ! reposition pointer to next row
               ENDIF
            ELSEIF(k==nz+1) THEN
               a(nNonZeros:nNonZeros)=half
               ja(nNonzeros:nNonZeros)= l-nx*ny
               ia(m)=nNonZeros
               nNonZeros=nNonZeros+1
            ELSE
               a(nNonZeros:nNonZeros+1)=(/half,half/)
               ja(nNonZeros:nNonZeros+1)=(/l-nx*ny,l/)
               ia(m)=nNonzeros
               nNonZeros=nNonZeros+1+1
            ENDIF

         ENDDO
      ENDDO
   ENDDO
nNonZeros=nNonZeros-1;ia(m+1)=nNonZeros+1

END SELECT

IF(size(a)/=nNonzeros) THEN
   PRINT*,'nNonzeros is not correct in e2c'
   STOP
ENDIF
IF(MAXVAL(ja)>nx*ny*nz) THEN
   PRINT*,'problems with definition of matrix in e2c'
   PRINT*,'columns exceed size'
   STOP
END IF

END SUBROUTINE c2e

SUBROUTINE gc2eb(a,ja,ia,nNonZeros,delta,dir,nx,ny,nz)
! Defines and fills in the sparse matrix corresponding to a gradient 
! in the direction specified by dir, applied to a cell-centered field 
! and returned on cell edges (c2e) along which the grid spacing is
! delta. The sparse matrix is (nx-I)*(ny-J)*(nz-K) X nx*ny*nz 
! (hence the b as in big)  where
! I=1 iff dir=1, J=1 iff dir=2 and K=1 iff dir=3 
! It returns the gradient of a field whose boundary conditions are known
! Declarations
INTEGER, INTENT (IN) :: nx,ny,nz,dir
REAL (KIND=Precision), INTENT (IN) :: delta
INTEGER, INTENT (OUT) :: nNonZeros
INTEGER, POINTER, DIMENSION (:) :: ia,ja
REAL (KIND=Precision), POINTER, DIMENSION (:) :: a

INTEGER :: i,j,k,l,m,d(3,3),stride(3),nrow
d(:,1)=(/1,0,0/);d(:,2)=(/0,1,0/);d(:,3)=(/0,0,1/)
stride(:)=(/1,nx ,(nx )*(ny )/)
! release memory
IF(ASSOCIATED(a)) THEN
   DEALLOCATE(a); NULLIFY(a)
ENDIF
IF(ASSOCIATED(ia)) THEN
   DEALLOCATE(ia); NULLIFY(ia)
ENDIF
IF(ASSOCIATED(ja)) THEN
   DEALLOCATE(ja); NULLIFY(ja)
ENDIF



nNonZeros=(nx -d(1,dir))*(ny -d(2,dir))*(nz -d(3,dir))*2
nrow=(nx -d(1,dir))*(ny -d(2,dir))*(nz -d(3,dir))
ALLOCATE(a(nNonzeros),ja(nNonzeros),ia(nrow+1))
nNonzeros=1
DO k=d(3,dir),nz-1!d(3,dir)
   DO j=d(2,dir),ny-1!+d(2,dir)
      DO i=d(1,dir),nx-1!+d(1,dir)
         l=(k)*(nx )*(ny )+(j)*(nx )+i+1 ! vectorized column index
         m=(k-d(3,dir))*(ny -d(2,dir))*(nx -d(1,dir))+&
              (j-d(2,dir))*(nx -d(1,dir))+(i-d(1,dir))+1
         a(nNonzeros:nNonzeros+1)=(/-1/delta,1/delta/)
         ja(nNonZeros:nNonzeros+1)=(/l-stride(dir),l/)
         ia(m)=nNonZeros
         nNonZeros=nNonZeros+1+1
      ENDDo
   ENDDO
ENDDO
nNonZeros=nNonZeros-1;ia(m+1)=nNonZeros+1

IF(size(a)/=nNonzeros) THEN
   PRINT*,'nNonzeros is not correct in gc2eb'
   STOP
ENDIF
IF(MAXVAL(ja)>(nx )*(ny )*(nz )) THEN
   PRINT*,'problems with definition of matrix in gc2eb'
   PRINT*,'columns exceed size'
   STOP
END IF


END SUBROUTINE gc2eb

SUBROUTINE c2eb(a,ja,ia,nNonZeros,dir,nx,ny,nz)
! Defines and fills in the sparse matrix corresponding to an average 
! in the direction specified by dir, applied to a cell-centered field 
! and returned on cell edges (c2e) along which the grid spacing is
! delta. The sparse matrix is (nx -I)*(ny -J)*(nz -K) X nx *ny *nz  
! (hence the b as in big)  where
! I=1 iff dir=1, J=1 iff dir=2 and K=1 iff dir=3 
! Note how the dimension along which the average is takes is reduced by one.
! Declarations
INTEGER, INTENT (IN) :: nx,ny,nz,dir
INTEGER, INTENT (OUT) :: nNonZeros
INTEGER, POINTER, DIMENSION (:) :: ia,ja
REAL (KIND=Precision), POINTER, DIMENSION (:) :: a

INTEGER :: i,j,k,l,m,d(3,3),stride(3),nrow
d(:,1)=(/1,0,0/);d(:,2)=(/0,1,0/);d(:,3)=(/0,0,1/)
stride(:)=(/1,nx ,(nx )*(ny )/)
! release memory
IF(ASSOCIATED(a)) THEN
   DEALLOCATE(a); NULLIFY(a)
ENDIF
IF(ASSOCIATED(ia)) THEN
   DEALLOCATE(ia); NULLIFY(ia)
ENDIF
IF(ASSOCIATED(ja)) THEN
   DEALLOCATE(ja); NULLIFY(ja)
ENDIF

nNonZeros=(nx -d(1,dir))*(ny -d(2,dir))*(nz -d(3,dir))*2
nrow=(nx -d(1,dir))*(ny -d(2,dir))*(nz -d(3,dir))
ALLOCATE(a(nNonzeros),ja(nNonzeros),ia(nrow+1))
nNonzeros=1
DO k=d(3,dir),nz-1!d(3,dir)
   DO j=d(2,dir),ny-1!+d(2,dir)
      DO i=d(1,dir),nx-1!+d(1,dir)
         l=(k)*(nx )*(ny )+(j)*(nx )+i+1 ! vectorized column index
         m=(k-d(3,dir))*(ny -d(2,dir))*(nx -d(1,dir))+&
              (j-d(2,dir))*(nx -d(1,dir))+(i-d(1,dir))+1
         a(nNonzeros:nNonzeros+1)=(/half,half/)
         ja(nNonZeros:nNonzeros+1)=(/l-stride(dir),l/)
         ia(m)=nNonZeros
         nNonZeros=nNonZeros+1+1
      ENDDo
   ENDDO
ENDDO
nNonZeros=nNonZeros-1;ia(m+1)=nNonZeros+1



IF(size(a)/=nNonzeros) THEN
   PRINT*,'nNonzeros is not correct in gc2eb'
   STOP
ENDIF
IF(MAXVAL(ja)>(nx )*(ny )*(nz )) THEN
   PRINT*,'problems with definition of matrix in gc2eb'
   PRINT*,'columns exceed size'
   STOP
END IF



END SUBROUTINE c2eb

SUBROUTINE strip(a,ja,ia,nNonZeros,dir,nx,ny,nz)
! Defines and fills in the sparse matrix corresponding to stripping 
! beginning and end points direction specified by dir, 
! applied to a cell-centered field. 
! The sparse matrix is (nx -2I)*(ny -2J)*(nz -2K) X nx *ny *nz  
! (hence the b as in big)  where
! I=1 iff dir=1, J=1 iff dir=2 and K=1 iff dir=3 

INTEGER, INTENT (IN) :: nx,ny,nz,dir
INTEGER, INTENT (OUT) :: nNonZeros
INTEGER, POINTER, DIMENSION (:) :: ia,ja
REAL (KIND=Precision), POINTER, DIMENSION (:) :: a
INTEGER :: i,j,k,l,m,d(3,3),stride(3),nrow
d(:,1)=(/1,0,0/);d(:,2)=(/0,1,0/);d(:,3)=(/0,0,1/)

IF(ASSOCIATED(a)) THEN
   DEALLOCATE(a); NULLIFY(a)
ENDIF
IF(ASSOCIATED(ia)) THEN
   DEALLOCATE(ia); NULLIFY(ia)
ENDIF
IF(ASSOCIATED(ja)) THEN
   DEALLOCATE(ja); NULLIFY(ja)
ENDIF

nNonZeros=(nx -2*d(1,dir))*(ny -2*d(2,dir))*(nz -2*d(3,dir))
nrow=(nx -2*d(1,dir))*(ny -2*d(2,dir))*(nz -2*d(3,dir))
ALLOCATE(a(nNonzeros),ja(nNonzeros),ia(nrow+1))
nNonzeros=1
DO k=d(3,dir),nz-1-d(3,dir)
   DO j=d(2,dir),ny-1-d(2,dir)
      DO i=d(1,dir),nx-1-d(1,dir)
         l=(k)*(nx )*(ny )+(j)*(nx )+i+1 ! vectorized column index
         m=(k-d(3,dir))*(ny -2*d(2,dir))*(nx -2*d(1,dir))+&
              (j-d(2,dir))*(nx -2*d(1,dir))+(i-d(1,dir))+1

         a(nNonzeros:nNonzeros)=(/1.d0/)
         ja(nNonZeros:nNonzeros)=(/l/)
         ia(m)=nNonZeros
         nNonZeros=nNonZeros+1
      ENDDo
   ENDDO
ENDDO
nNonZeros=nNonZeros-1;ia(m+1)=nNonZeros+1



IF(size(a)/=nNonzeros) THEN
   PRINT*,'nNonzeros is not correct in gc2eb'
   STOP
ENDIF
IF(MAXVAL(ja)>(nx )*(ny )*(nz )) THEN
   PRINT*,'problems with definition of matrix in gc2eb'
   PRINT*,'columns exceed size'
   STOP
END IF

END SUBROUTINE strip

SUBROUTINE zeroes(a,ja,ia,nNonZeros,dir,nx,ny,nz)
! Defines and fills in the sparse matrix corresponding to zeroing
! the edges 
! beginning and end points direction specified by dir, 
! applied to a cell-centered field. 
! The sparse matrix is (nx )*(ny )*(nz ) X nx *ny *nz  
! (hence the b as in big) since we actually keep the zeros

INTEGER, INTENT (IN) :: nx,ny,nz,dir
INTEGER, INTENT (OUT) :: nNonZeros
INTEGER, POINTER, DIMENSION (:) :: ia,ja
REAL (KIND=Precision), POINTER, DIMENSION (:) :: a
INTEGER :: i,j,k,l,m,d(3,3),stride(3),nrow
d(:,1)=(/1,0,0/);d(:,2)=(/0,1,0/);d(:,3)=(/0,0,1/)

IF(ASSOCIATED(a)) THEN
   DEALLOCATE(a); NULLIFY(a)
ENDIF
IF(ASSOCIATED(ia)) THEN
   DEALLOCATE(ia); NULLIFY(ia)
ENDIF
IF(ASSOCIATED(ja)) THEN
   DEALLOCATE(ja); NULLIFY(ja)
ENDIF

nNonZeros=nx*ny*nz
nrow=nNonZeros
ALLOCATE(a(nNonzeros),ja(nNonzeros),ia(nrow+1))
a=1.d0

SELECT CASE (dir)
CASE (1)
DO k=1,nz
   DO j=1,ny
      DO i=1,nx
         l=(k-1)*(nx )*(ny )+(j-1)*(nx )+i ! vectorized column index
         IF(i==1.OR.i==nx) a(l)=zero
         ja(l)=l
         ia(l)=l
      ENDDo
   ENDDO
ENDDO
CASE (2)
DO k=1,nz
   DO j=1,ny
      DO i=1,nx
         l=(k-1)*(nx )*(ny-1 )+(j)*(nx )+i ! vectorized column index
         IF(j==1.OR.j==ny) a(l)=zero
         ja(l)=l
         ia(l)=l
            
      ENDDo
   ENDDO
ENDDO
CASE (3)
DO k=1,nz
   DO j=1,ny
      DO i=1,nx
         l=(k-1)*(nx )*(ny )+(j-1)*(nx )+i ! vectorized column index
         IF(k==1.OR.k==nz) a(l)=zero
         ja(l)=l
         ia(l)=l
            
      ENDDo
   ENDDO
ENDDO

END SELECT

ia(l+1)=l+1;
IF(size(a)/=nNonzeros) THEN
   PRINT*,'nNonzeros is not correct in gc2eb'
   STOP
ENDIF
IF(MAXVAL(ja)>(nx )*(ny )*(nz )) THEN
   PRINT*,'problems with definition of matrix in gc2eb'
   PRINT*,'columns exceed size'
   STOP
END IF

END SUBROUTINE zeroes

SUBROUTINE zeroi(a,ja,ia,nNonZeros,dir,nx,ny,nz)
! Defines and fills in the sparse matrix corresponding to a zero operator 
! in the direction specified by dir, applied to a cell-centered field 
! and returned on cell edges (c2e) along which the grid spacing is
! delta. The sparse matrix is (nx -I)*(ny -J)*(nz -K) X nx *ny *nz  
! (hence the b as in big)  where
! I=1 iff dir=1, J=1 iff dir=2 and K=1 iff dir=3 
! Note how the dimension along which the average is takes is reduced by one.
! Declarations
INTEGER, INTENT (IN) :: nx,ny,nz,dir
INTEGER, INTENT (OUT) :: nNonZeros
INTEGER, POINTER, DIMENSION (:) :: ia,ja
REAL (KIND=Precision), POINTER, DIMENSION (:) :: a

INTEGER :: i,j,k,l,m,d(3,3),stride(3),nrow
d(:,1)=(/1,0,0/);d(:,2)=(/0,1,0/);d(:,3)=(/0,0,1/)
stride(:)=(/1,nx ,(nx )*(ny )/)
! release memory
IF(ASSOCIATED(a)) THEN
   DEALLOCATE(a); NULLIFY(a)
ENDIF
IF(ASSOCIATED(ia)) THEN
   DEALLOCATE(ia); NULLIFY(ia)
ENDIF
IF(ASSOCIATED(ja)) THEN
   DEALLOCATE(ja); NULLIFY(ja)
ENDIF

nNonZeros=(nx -d(1,dir))*(ny -d(2,dir))*(nz -d(3,dir))
nrow=(nx -d(1,dir))*(ny -d(2,dir))*(nz -d(3,dir))
ALLOCATE(a(nNonzeros),ja(nNonzeros),ia(nrow+1))
nNonzeros=1
DO k=d(3,dir),nz-1!d(3,dir)
   DO j=d(2,dir),ny-1!+d(2,dir)
      DO i=d(1,dir),nx-1!+d(1,dir)
         l=(k)*(nx )*(ny )+(j)*(nx )+i+1 ! vectorized column index
         m=(k-d(3,dir))*(ny -d(2,dir))*(nx -d(1,dir))+&
              (j-d(2,dir))*(nx -d(1,dir))+(i-d(1,dir))+1
         a(nNonzeros:nNonzeros)=(/zero/)
         ja(nNonZeros:nNonzeros)=(/l/)
         ia(m)=nNonZeros
         nNonZeros=nNonZeros+1
      ENDDo
   ENDDO
ENDDO
nNonZeros=nNonZeros-1;ia(m+1)=nNonZeros+1



IF(size(a)/=nNonzeros) THEN
   PRINT*,'nNonzeros is not correct in gc2eb'
   STOP
ENDIF
IF(MAXVAL(ja)>(nx )*(ny )*(nz )) THEN
   PRINT*,'problems with definition of matrix in gc2eb'
   PRINT*,'columns exceed size'
   STOP
END IF



END SUBROUTINE zeroi

SUBROUTINE dress(a,ja,ia,nNonZeros,dir,nx,ny,nz)

! Defines and fills in the sparse matrix corresponding to the 
! operator that extends a field to include ghost points (set to zero)  
! in the direction specified by dir, applied to a cell-centered field 
! and returned on cell edges (c2e) along which the grid spacing is
! delta. The sparse matrix is (nx +2*)*(ny +2J)*(nz +2K) X nx *ny *nz  
! (hence the b as in big)  where
! I=1 iff dir=1, J=1 iff dir=2 and K=1 iff dir=3. If dir=0 then I=1,J=1,K=1 

! Declarations
INTEGER, INTENT (IN) :: nx,ny,nz,dir
INTEGER, INTENT (OUT) :: nNonZeros
INTEGER, POINTER, DIMENSION (:) :: ia,ja
REAL (KIND=Precision), POINTER, DIMENSION (:) :: a


INTEGER :: i,j,k,l,m,d(3,0:3),stride(3),nrow
d(:,1)=(/1,0,0/);d(:,2)=(/0,1,0/);d(:,3)=(/0,0,1/)
d(:,0)=(/1,1,1/)


! release memory
IF(ASSOCIATED(a)) THEN
   DEALLOCATE(a); NULLIFY(a)
ENDIF
IF(ASSOCIATED(ia)) THEN
   DEALLOCATE(ia); NULLIFY(ia)
ENDIF
IF(ASSOCIATED(ja)) THEN
   DEALLOCATE(ja); NULLIFY(ja)
ENDIF

nrow=(nx+2*d(1,dir))*(ny+2*d(1,dir))*(nz+2*d(3,dir))
nNonZeros=nx*ny*nz
ALLOCATE(a(nNonzeros),ja(nNonzeros),ia(nrow+1))
nNonzeros=1
DO k=1-d(3,dir),nz+d(3,dir)
   DO j=1-d(2,dir),ny+d(2,dir)
      DO i=1-d(1,dir),nx+d(1,dir)
         l=(k-1)*nx*ny+(j-1)*nx+i ! vectorized column index
         m=(k-(1-d(3,dir)))*(ny+2*d(2,dir))*(nx+2*d(1,dir))+(j-(1-d(2,dir)))*(nx+2*d(1,dir))+(i-(1-d(1,dir)))+1 ! vectorized row index  
         ia(m)=nNonZeros
         
         IF(MOD(k,nz+1)*MOD(j,ny+1)*MOD(i,nx+1)/=0) THEN ! interior points

            a(nNonZeros:nNonZeros)=(/one/)
            ja(nNonZeros:nNonZeros)=(/l/)
            
            nNonZeros=nNonZeros+1
         ENDIF
      ENDDO
   ENDDO
ENDDO
nNonZeros=nNonZeros-1;ia(m+1)=nNonZeros+1
IF(size(a)/=nNonzeros) THEN
   PRINT*,'nNonzeros is not correct in gc2eb'
   STOP
ENDIF
IF(MAXVAL(ja)>(nx )*(ny )*(nz )) THEN
   PRINT*,'problems with definition of matrix in gc2eb'
   PRINT*,'columns exceed size'
   STOP
END IF
END SUBROUTINE dress

!***********************
END MODULE Basic_matrix
!***********************
